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Abstract. It is well known that multigrid methods are optimally efficient for 
solution of elliptic equations (O(TV)), which means that effort is proportional 
to the number of points at which the solution is evaluated). Thus this is an 
ideal method to solve the initial data/constraint equations in General Relativity 
for (for instance) black hole interactions, or for other strong-field gravitational 
configurations. Recent efforts have produced finite difference multigrid solvers 
for domains with holes (excised regions). We present here the extension of 
these concepts to higher order (fourth-, sixth- and eigth-order). The high order 
convergence allows rapid solution on relatively small computational grids. Also, 
general relativity evolution codes are moving to typically fourth-order; data have 
to be computed at least as accurately as this same order for straightfoward 
demonstration of the proper order of convergence in the evolution. 

Our vertex-centered multigrid code demonstrates globally high-order-accurate 
solutions of elliptic equations over domains containing holes, in two spatial 
dimensions with fixed (Dirichlet) outer boundary conditions, and in three spatial 
dimensions with Robin outer boundary conditions. We demonstrate a "real world" 
3-dimensional problem which is the solution of the conformally flat Hamiltonian 
constraint of General Relativity. The success of this method depends on: a) the 
choice of the discretization near the holes; b) the definition of the location of the 
inner boundary, which allows resolution of the hole even on the coarsest grids; 
and on maintaining the same order of convergence at the boundaries as in the 
interior of the computational domain. 



1. Introduction 

Solving Einstein's equation in 3+1 form [TJ [5] requires that a set of elliptic (or 
quasi-elliptic) constraint equations be satisfied at all times. The Bianchi identities 
ensure that, given a set of initial data which analytically satisfy the constraints, 
the subsequent analytically evolved variables will also satisfy the constraints. In 
numerical solutions of Einstein's equations, however, the constraints are not preserved 
exactly. Thus errors will arise as the simulation proceeds, and the extent to which the 
numerical solutions actually reflect the true solutions of the analytic equations is an 
ongoing area of research (e.g., [H SJ [3 [7] ) . Apart from the issue of the accuracy 
of the solutions obtained, there is also the problem of numerical (in)stability, may in 
some cases be related to the lack of preservation of the constraints. There is active 
interest in constrained general relativity evolution schemes (which repeatedly solves 
the constraint equations during the evolution [51 151 HP] ) . 

The multigrid method is a particularly attractive method to solve elliptic 
equations because it is an optimal method, i.e. it requires only O(N) operations, 
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where N is the number of unknowns (proportional to the number of locations at 
which knowledge of the independent variable is required) . Furthermore, the multigrid 
method is not especially difficult to implement. We study the multigrid problem on 
domains with holes, because while many formulations now can produce moving black 
holes without excising the holes, evolution with excision in a necessary requirement 
for fully studying and understanding black hole evolution. For instance, posing data 
for black holes with Kerr parameter a very near unity are straightforward in excised 
data, but very difficult to achieve in non-excised data (such as puncture data[TT]). 

It has long been thought infeasible in practice to obtain generic results which 
have everywhere the same order of accuracy as the finite difference scheme employed, 
for domains containing holes. This argument was refuted in [12] for 3D solutions of 
second-order accuracy. This paper addresses the problem for domains with coordinate- 
spherical holes in the context of higher order accuracy. Multigrid provides a simple, 
robust method to achieve solutions in such cases. 

If the multigrid solver is to be used to provide initial data for a numerical evolution 
code, the error of the initial solution need only be below the truncation error of the 
finite difference scheme used for evolution [13]. However high order methods are still 
extremely valuable in this context because they require much smaller computational 
resources for the same accuracy. Thus this high order approach may dramatically 
lower the elliptic - solution load in constrained evolutions. Furthermore, to correctly 
demonstrate the convergence of evolution codes (now typically written in terms of 
fourth-order stencils), one must have data of at least the same order of convergence. 

Finite difference techniques to solve the constraint equations are a standard 
activity in relativity, and the properties of the constraint equations are well known 
(for the Euclidean, constant-mean-curvature background and for the signs of the 
solutions we use here QTj [15j [16]). We focus on the benefits of high-order multigrid 
implementation. The validity of our results will be established by demonstrating 
controlled high order convergence to known exact solutions. 

Section [2] assumes familiarity with the multigrid method, but a small overview is 
presented. Section [3] describes our new scheme for inner boundary points and presents 
the results for a simple problem in a 2D domain. Section 2] extends this to 3D. Section 
[5] summarizes our results. 

2. Overview of Multigrid 

Application of the multigrid method [17j . on domains with holes — or on domains 
with "irregular boundaries" in general — has received only modest attention [18] , [19] ; 
these are cell centered finite difference codes. The "BAM" code [20] provides access 
to the multigrid method in numerical relativity. It is a vertex centered second-order 
code and features the ability to handle domains with multiple holes, though with 
restrictions on hole placement. Hawley and Matzner [T2] introduced a second-order 
vertex-centered multigrid code that accommodates arbitrary placement of excised 
holes. The innovation in this current paper is solution to higher convergence order (we 
give fourth- and sixth- order examples, and some preliminary, eighth-order results) 
which dramatically reduces the time and storage requirements of the elliptic solution. 

The multigrid scheme [17] has received considerable attention in the literature, 
and is the subject of numerous articles, conferences, reviews and books (e.g., 
[211 [22l [23l [24l [25l [26]). It is essentially a clever means of eliminating successive 
wavelength-components of the error via the use of relaxation at multiple spatial scales. 
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Here we give a very brief overview of the multigrid method, following the notes by 
Choptuik [27] . (Introductions to multigrid applications in numerical relativity are also 
found in Choptuik and Unruh [25J and Brandt [25].) We want to solve a continuum 
differential equation Cu = /, where £ is a differential operator, / is some right hand 
side, and u is the solution we wish to obtain. We discretize this differential equation 
into a difference equation on some grid (or lattice) with uniform spacing h: 

&u h = f h , (1) 

where u h is the exact solution of this discrete equation, and lim^o u — u - The 
discretization h refers to the finest grid, of a hierarchy of vertex-centered grids. Each 
grid at multigrid level I is a square lattice having 2 l + 1 grid points along each edge. 
The grids have uniform spacing hi = 2~ l in both x and y directions, and the grid 
points are denoted with indices i and j in the x and y directions, respectively, e.g., 
u(ihi,jhi) ~ u\y 

Rather than attempting to solve Eq. (fT]) directly via a costly matrix inversion, 
we apply an iterative solution method. At any step in the iteration, we have an 
approximate solution u h ~ u h . In this iterative algorithm, the original guess u^ ld is 
brought closer to u h by applying some (approximate) correction: 

u^:=^ ld + v h (2) 

For nonlinear operators the correction is via the Full Approximation Storage 
(FAS) method [H EES E3J [2J [2J [2J (2J [2J . 



2.1. V- Cycles and the Full Multigrid Algorithm 

The solution algorithm takes the form of a V-cycle, in which we start with an initial 
guess on the fine grid, at multigrid level l max (the finest grid has 2' max + 1 vertices 
along each edge). Then we perform some number of smoothing sweeps. A smoothing 
sweep is one iteration of a relaxation solver for the elliptic equation. The effect of 
such a step is to reduce the short wavelength error on the grid, i.e. it "smooths" the 
approximate solution. At this point we wish to update Eq.([2]). 

To accomplish this the multigrid method introduces I? , the restriction operator, 
which maps values from the fine grid to the next coarser grid via some weighted 
averaging operation, and , an interpolation or prolongation operator, which maps 
values from a coarse grid to the next finer grid via some interpolation operation. In 
the Full Approximation Storage multigrid method we express the correction as 

ul w :=u% M + lZ h (u 2h -I 2 h h u h ). (3) 

Eq([3]) defines a coarse grid correction (CGC). 

In Eq([3]), u 2h is the exact solution to {1} on the next coarser grid. However, 
we may approximately solve for the coarser grid u 2h ~ u 2h by repeating Eq. ([3]) 
on the coarse grid, which then refers to the next (even) coarser grid. We continue 
smoothing and restricting to coarser grids until we arrive at a grid coarse enough 
to solve the resulting coarse grid equation 'exactly' (i.e., to machine precision), at 
minimal computational cost. Equation ([3]) can then be used to correct the next 
finer grid solution. We can thus iteratively correct the solutions on each next finer 
grid, perhaps with additional smoothing operations performed before moving to each 
finer grid. One may carry out a number of smoothing sweeps at a given refinement 
level before proceeding to solve on the next coarser grid (pre-sweeps), and a (perhaps 
different) number of post-sweeps after solving on the next coarser grid. We use the 
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same numbers of pre- and the same number of post- sweeps at every level, though the 
numbers could in principle be varied at different levels. 

On all grids except the coarsest grid, we only smooth the error, and we solve the 
difference equation exactly only on the coarsest grid. In practice we solve the coarse 
grid difference equation by relaxation, which is cheap on the coarsest grid. The entire 
process, as described is called a V-cycle. A full solution may consist of one or more 
V- cycles. 

3. Solution of a Nonlinear Poisson Equation in 2D 

We solve the equation 

d 2 d 2 

— u(x, y) + q^u{x, y) + <Ju 2 (x, y) = f(x, y), (4) 

on a domain £1 with coordinate ranges [0, 0] to [1, 1], and subject to Dirichlet conditions 
at the outer boundary: u(x,y)\gn = 0. The function f(x,y) is chosen such that the 
solution is 

u{x, y) = sin(7ra;) sin(7n/), (5) 

Inner Boundary Conditions 

We have added features to handle holes in the domain. In [TJ] the inner boundary 
is given by the outer edge of the points comprising the excision mask. This has the 
consequence that the "size" of the excision region (i.e., the area of the convex hull of 
the data points comprising the mask) on finer grids is always equal to or greater than 
the size of the excision region on coarser grids. Here we present a novel modification 
of the inner boundary algorithm: at the inner boundary we insert a new set of points 
into the grid, and apply Dirichlet boundary conditions (the exact solution) to these 
points. The boundary conditions are applied at each level. This has the effect of 
keeping the size of the excised region essentially constant at all levels. 

These extra points are added at each level at every location at which a grid line 
intersects the exact inner boundary location. The new boundary points specify the 
exact solution to the inner boundary conditions. The neighboring grid points use this 
new point in their smoothing and stencil operations taking care to account for the 
irregular grid size caused by adding this point. Figure Q] illustrates this. Note that 
this has the important result that the hole is resolved on even the coarsest grid. We 
will find that this approach dramatically reduces the solution error in domains with 
holes. Although Figure [T] shows a circular excision centered on a grid point, we have 
verified that convergence is maintained even for arbitrary placement of the excision 
circle (or sphere, in 3D). 

Smoothing Operations 

Except on the coarsest grid, we do not attempt to solve the difference equation 
Instead we use a few sweeps of a relaxation method. This is called smoothing. A "red- 
black" Gauss-Seidel Newton iteration is used to relax the solution(smooth the error). 
For all points that have a regular spacing h with their neighbors (interior points) this 
is the traditional method. The values of h differ by a factor of 2 between "adjacent" 
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Figure 1. Extra grid points points added at intersection of the inner boundary 
and the grid lines. Dirichlet boundary conditions are used to fill in value for these 
points. 



grid levels. For these 2D examples we use a second-order discretization. On interior 
points, Eq([4]) is discretized as 



r GS = h + u%-i,j + u id+1 + u itj -i - ^Uij) 



™h - hi ( 6 ) 



-now _ -old r GS 

« 2a,,,,, \l, () 

But for points where the spacings are not regular, we introduce new notation here. 

Define h i+ 1 ■ as the spacing in the x direction between points at i,j and i + 

similarly h^i j is the spacing between i,j and i — Note that there can be 

multiple substitutions for an excised grid point, so we choose a particular one, and 

the associated h j+ i 
2 >j 



ras 
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<?«:.,, f.., (8) 

1 1 , 2 



+ { T^ + T^^1 h < 9 > 

Note that because of the unequal spacing, Eqs([8l- [T0)) are in fact only first-order near 
the inner excision boundary. For 2D we find that this does not contaminate the high 
order convergence of the solutions. 

Restriction and Prolongation Operators 

The restriction operator I^ h we use is the so-called "half-weighted" average on 
normal interior points, in which coarse grid values (indexed by / and J for clarity) are 
a weighted average of the fine grid values over a nearby region of the physical domain: 

u] h j = IT u h = \u h id + 1 + ul ltj + u h i>j+1 + u h id _ x ] , (11) 

where i = 21 —1 and j = 2 J — 1. 

For the prolongation operator ! we use simple bilinear interpolation. 

3.1. 2D Results 



We solve EqQ with a = 1. We take l max = 7 (the finest grid has 2 7 + 1 = 129 
vertices along an edge), and l m i n = 2 (the coarsest grid has 2 2 + 1 = 5 vertices along 
an edge) . Figure [2] is a plot of the error over the entire domain, with no holes in the 
domain. Most implementations of multigrid with irregular inner boundaries end up 
with higher errors than ones without the holes, but the extra inner boundary condition 
might also help "tie the solution down" . Figure [3] shows solution error for a domain 
with an excised region (and inner boundary) of radius r = 0.129. The scales of Figures 
[2] through [4] are matched, so it is easy to see that the error decreases. Also note that 
the error is very tightly bounded close to the inner boundary; evidently the inner 
boundary condition is reducing the maximum error in the domain. Figure |4] which 
shows the solution error when the domain contains two holes, confirms this. 

Figure [5] shows the average error over the domain for 1 and 2 V-cycles (2 pre 
and 2 post smoothing runs) with and without a hole. We need 2 V-cycles to achieve 
second-order convergence. As noted above, the error for the runs with holes is smaller 
than in the the base runs (runs without holes). 

Figure [6] is a similar plot comparing the base runs with runs over domains with 
two holes. Again the error is lower than the base run without holes. 

Figures [2] through [6] demonstrate the dramatic improvement obtained by our new 
inner boundary treatment. In particular Figure [3] describes the error for precisely 
the same problem as that treated in Figure 3 of ref [12] . The only difference is the 
inner boundary treatment. In (,12J, Figure 3) the error fluctuates and the maximum 
occurs near the hole; maximum error ~ 0.7 x 10~ 3 . The error near the inner boundary 
also shows noticeable imprinting from the cartesian grid. This is achieved for eight 
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Figure 2. A plot of the solution error e = \u — u h \ for the 2D case. We present 
a second-order solution with no holes to provide a metric for comparison for the 
rest of the results. i ma x = 7 (the finest grid has 2 7 + 1 = 129 vertices along an 
edge), and Z m i n = 2 (the coarsest grid has 2 2 + 1 = 5 vertices along an edge). 2 
V-cycle and 2 pre- and 2 post-CGC smoothing sweeps. 

multigrid levels (finest grid of 2 s + 1 = 257 vertices along an edge). In this present 
work, Figure ([3]), the error near the hole is quite small and smooth and shows no 
noticeable rectangular character. The maximum error across the grid is ~ 0.1 x 10~ 4 
(more than an order of magnitude smaller than that of [II]!). This is accomplished 
with only seven multigrid levels. 

Clearly the new inner boundary treatment (inserting new points), with the 
important feature of resolving the hole on every grid, dramatically outperforms 
the previous inner boundary method, which applied boundary conditions only at 
regular grid points (thus giving a "rougher" - more "lego-like" definition of the inner 
boundary). 

While there are many ways this work may be improved, it is encouraging that 
even without mesh refinement it represents a dramatic improvement in accuracy. 

4. 3D Simulations 

The Hamiltonian constraint for a single black hole in a conformally flat background 
geometry yields the following equation: 

A 2 

V 2 u(x, y, z) - K 2 u b {x, y, z) + — , = f(x, y, z), (12) 

u'(x,y,z) 

This is Poisson equation with two nonlinear terms: V 2 is the usual Laplacian in 
Euclidean space. K 2 and A 2 are arbitrary positive real constants related to the rate 
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of expansion for the 3-space for which (|12p is the constraint equation, fix, y, z) is 
related to the energy density in the 3-space, and can be chosen such that the resulting 
u{x, y, z) has some known (exact) form by which we can check our numerical results. 
We adjust the matter source term f(x,y,z) such that the solution is 1 + 2M/r. We 
choose M = 1 and K, A = 1. 

For all the 3D simulations we use the domain [—5.0, —5.0, —5.0] to [5.0, 5.0, 5.0] 
and a hole of radius r max = 1.29. 

For the outer boundary in the 3D code to solve the conformally flat background, 
we apply the Robin condition. The Robin boundary condition is 

djrju - 1)) = Q 
dr 

This Robin condition is standard condition used in relativity ([3~2])- We choose to 
follow Alcubierre [30] and [12] and take derivatives only in direction normal to the 
faces of our cubical domain. Since the stencil for the differencing about a boundary 
point is not symmetric about the boundary point we need an extra point to maintain 
second order accuracy. We do so for second order differencing, and for higher order 
runs we include the required number points in the stencil to match the order to that 
of the interior points. 



Restriction and Prolongation Operators 
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Figure 5. Average error with and without a central hole in the 2D case. The 
two "no hole" lines lie atop one another. 
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Figure 6. Plot of average error with and without two holes in the 2D case. The 
two "no hole" lines lie atop one another. 



We use the half-weighted restriction operator on normal interior points, with a 
weight of 1/12 for each neighborhood point along each of the x, y, z direction. 

"?fj,Jf = iff 1 U h = 2™i,j,k + {^i+l,j,k + 'Ui-l,j,k Jr 

"i,3+l,* + "ij-l.fc + ^i,j,k+l + ^i,j,k-l) ' U^) 

where i = 21 — 1 and j = 2J — 1. For points near the boundary, we adjust the weight 
to take account of the different spacing. 

More complicated restriction operators are oft en ( 12\) used which include 
weighing all the points on a cube centered around the coarse grid point. The results of 
our simulation are not sensitive to the restriction operator used, so we prefer the simple 
half-weighted one described above. For the prolongation operator l£ h , we continue to 
use trilinear interpolation. 

4-1. Second- order convergence 



Smoothing Operator 

Again similar to the 2D case, we have to make exceptions for the interior points near 
the boundary. With regular spacing the usual estimate used to calculate V 2 is done 
as follows; thus in each direction, (say x), the second derivative is 

—J = h~ 2 (u l+hj + u^ hj - 2u it j) + 0(h 2 ). (14) 
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For irregular spacings a simple modification would be 
d 2 u 2 



0(h). (15) 



However because of the irregular spacing in Eq. (|15p the error during smoothing and 
hence the final solution near the inner boundary will be first-order instead of second. 
(This problem should not exist for traditional implementations that maintain a regular 
grid spacing.) Unlike the 2D case, we have found it essential in the 3D case to have 
consistent differencing order at the boundary; otherwise smooth convergence at the 
desired order could not be achieved. Thus, in order to maintain second-order error 
bounds near the boundary, we include an extra interior point. So a total of four points 

(u(x),u(x + h),u(x + 2h),u(x + 3/i)) are used to calculate the second derivative in 

ft 2 ~ 

each of x, y, z coordinates. Then we solve the Taylor expansions for 

6(x+Al) = i(x) + g ftl + gM + g| +0W) (16) 
«c+w-«w^+S^St +0 «' (17) 

»(x + h 3 ) = u(x) + + _ J + _.JL + 0{ht) (18) 
As we can see in fig [7] we maintain second-order convergence on all points. 



4-2. Fourth-order convergence 

Given the success in second-order convergence we now attempt to extend this to 
fourth order. We find successful fourth-order convergence by increasing the size of the 
smoothing stencil to five points in each [x,y, z) direction: 

d 2 u h~ 2 

-Q^a = ~jy ( 16 ("'+ij + "*-i,3') ~ + Ui-2,j) - 30 "ij) + 0(h A ). (19) 

Notice that we increase the order of accuracy of the smoothing operator only. We do 
not modify the interpolation or prolongation operators. 

For points near the boundary we need to include one additional point, and 
describing the Taylor expansions. The Taylor equations for the values over a total 
of six points can then be solved for 0^ to 0(h 4 ) accuracy. 

The use of stencils that include these additional points could create problems with 
the the red-black sweeping of the Gauss smoothing. Points updated in the red sweep 
can immediately be used in calculations for other points within the same sweep. Thus 
with these additional points Gauss smoothing can introduce a bias in the direction of 
the sweep that could affect the error. We in fact observed this bias. We correct this 
bias as follows, in a sort of red-black Jacobi method. We assume a red-black labeling 
of the vertices (every nearest neighbor of a red is a black, every nearest neighbor of a 
black is a red) . We then compute the update of each red vertex, using as many points 
as appropriate to the the discretization order, but these updated red values are stored 
in a separate auxiliary grid (rather than being written back into the original grid as 
would be done in the Gauss procedure). After all red points are updated, they are 
all copied back into the original grid. Then the black points are similarly updated, 
stored in an auxiliary grid, and then copied back into the original grid. The extra 
memory required to hold intermediate results is not a concern since other steps in the 
multigrid always require extra work arrays, which are available for use at this time. 
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Notice that we increase the order of accuracy of the smoothing operator only. 
We do not modify the interpolation or prolongation operators. A further complication 
is that with larger stencils of high-order methods, there may not be enough "room" 
on the coarsest grid to construct the stencil. An obvious solution is to ensure that 
the coarsest grid is sufficiently fine that that the stencils do fit. Instead, we use the 
following approach: If there is insufficient room to construct a stencil at the desired 
order, we drop back to the next highest order that does fit. In practice this means the 
solution on the coarsest grid is always second (or even first) order. 

Regardless of the simplifications, increasing only the order of the smoothing, 
and even dropping back to low-order differencing on coarse grids, this approach has 
yielded fourth-order solutions. Figure [8] shows the fourth-order convergence we get 
with the combination of the new stencils and modified red-black smoothing. Figure [9] 
shows a comparison between the second-order and fourth-order methods of the error 
norms produced. We do need a larger number of V-cycles and/or smoothing cycles 
to get fourth-order convergence, compared to the second-order case. We have not 
yet attempted to reduce or optimize this, but with fourth-order convergence the error 
behavior is so much better than with lower-order schemes that the extra time spent 
on additional V-cycles is not a matter of concern. 

A plot of the errors across the entire domain is useful to visualize the performance 
of the second-order and fourth-order methods. Figures ITD1 and [TT1 are plot of the error 
across the central slice of the domain. The error is highest close to the inner boundary 
but the fourth-order runs contain the error to a very small neighborhood of the inner 
boundary, and the maximum error is a factor ~ 10 2 smaller than that for the second- 
order solution. 

4-3. Sixth-, and higher- order convergence 

The work was further extended to convergence of higher order than fourth. We include 
additional points into the stencil to explore sixth- and eighth-order convergence. The 
results are shown in Figure [131 We see that, if a sufficiently large number of levels 
is included, good convergence behavior is obtained to eighth order. It appears that 
the number of refinement levels to achieve convergence increases with the convergence 
order. 

5. Conclusions and further work 

The main ideas presented in this paper are 

• We augment the grid for each level by 

— Finding the intersection between every grid line at each level with the inner 
boundaries (holes). 

— At each of these points we insert a additional boundary grid point, whose 
solution in the Dirichlet boundary conditions is known a priori. 

We have demonstrated explicitly in the second-order 2D implementation that 
this procedure leads to much improved error at the inner boundary. The essential 
feature appears to be resolving the hole on every grid. 

• Since the grid spacings near the inner boundaries are now not fixed, we modified 
the Gauss-Seidel Newton iteration (used to smooth the error) to take one 
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Figure 7. Error norms. Stencils are chosen for second-order convergence and 
plots with different numbers of V-cycles, pre and post smoothings are presented. 
The two 4 Vcycle average error plots lie atop one another. 




v 1 Vcycle, 2 pre, 2 post 

1 Vcycle, 4 pre, 4 post 
+ 2 Vcycle, 2 pre, 2 post 

2 Vcycle, 4 pre, 4 post 
• 4 Vcycle, 2 pre, 2 post 
■ 4 Vcycle, 4 pre, 4 post 

□ □ 4 Vcycle, 4 pre, 4 post, max error 

3 4 5 



Figure 8. Error norms. The stencils chosen here should give fourth-order 
convergence. 
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Figure 9. Comparison for second-order and fourth-order, of average error, and 
the maximum error over the entire domain, which is always located near the inner 
boundary. The plot clearly shows that error at the boundary is also controlled to 
the correct order 



additional neighboring grid point in each of x, y, z coordinates to maintain the 
relaxation accuracy. 

• We achieve good convergence for a second-order implementation. And, by 
increasing the order of the smoothing operator (only) and even with a drop back to 
lower order on the coarsest grids, we find good convergence at fourth- and sixth- 
order. We have demonstrated a remarkably simple implementation to provide 
fourth- and sixth-order elliptic solutions. 

Problems with this approach are 

• One of the planned applications of the code is in constrained evolution. The 
manner in which we introduce grid points leads to some grid spacings near the 
boundary that are much smaller than the regular grid spacing. This may cause 
violations of the Courant-Friedrichs-Levy(CFL) constraint on the time steps for 
evolution. 

• The complexity of the data structures and of the coding increases for higher order 
so that parallelizing this code without introducing bugs might be difficult. 

Some of the results of this work have not yet been investigated to completely 
understand the behavior. For example, we do not fully understand why we need 4 
V-cycles for convergence in the 3D runs, nor have we fully understood the number of 
levels required to insure higher-order convergence. 
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Figure 11. Plot of error over a central slice x = 0.0. Fourth- order run. l D 



: 2, 4 V-cycle and 4 pre- and 4 post-CGC smoothing sweeps. 
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Figure 12. Plot of error norms. 
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Figure 13. Comparison of average error, and the maximum error over the entire 
domain, which is always located near the inner boundary, the plot clearly shows 
that error at the boundary is also controlled to consistent order. 
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